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ABSTRACT 

We use a model of polarized Galactic emission developed by the the Planck collab- 
oration to assess the impact of foregrounds on B-mode detection at low multipoles. 
Our main interest is to applications of noisy polarization data and in particular to 
u: assessing the feasibility of -B-mode detection by Planck. This limits the complexity of 

\ foreground subtraction techniques that can be applied to the data. We analyze inter- 

O |. nal linear combination techniques and show that the offset caused by the dominant 

q _E-mode polarization pattern leads to a fundamental limit of r ~ 0.1 for the tensor- 

s' . scalar ratio even in the absence of instrumental noise. We devise a simple, robust, 

' template fitting technique using multi- frequency polarization maps. We show that 

C$ . template fitting using Planck data alone offers a feasible way of recovering primordial 

£?-modes from dominant foreground contamination, even in the presence of noise on 
the data and templates. We implement and test a pixel-based scheme for computing 
^ ' the likelihood function of cosmological parameters at low multipoles that incorporates 

. foreground subtraction of noisy data. 

Key words: Methods: data analysis, statistical; Cosmology: cosmic microwave back- 
^J. , ground, large-scale structure of Universe 
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INTRODUCTION 



In the last decade, observations of cosmic microwave background (CMB) anisotropies have provided one of the most pow- 
erful probes of cosmology. By combining CMB anisotropy data with a variety of other data, many of the key parameters 
that define our Universe have been determined with unprecedented precision (see e.g. Komatsu el ai, 2008, and references 
' therein). Nevertheless, many important questions remained unanswered. One of the most important is the amplitude of a 
'B-mode' polarization signature in the CMB. Scalar perturbations generate during inflation generate purely a divergence-like 
B-mode polarization pattern in the CMB, whereas tensor perturbations would produce a distinctive curl-like (B-mode) polar- 
ization signature together with an B-mode pattern of roughly equal amplitude*(Zaldarriaga and Seljak, 1997; Kamionkowski, 
Kosowsky and Stebbins, 1997). 

A detection of a primordial B-mode anisotropy would provide crucial evidence that inflation actually took place. Fur- 
thermore, a measurement of the relative amplitude of the tensor and scalar primordial power spectra (the tensor-scalar ratio 
r, see Peiris et ai, 2003, for a precise definition) would fix the energy scale of inflation via 

U 1/4 » 3.3 x 10 16 r 1/4 GeV, (1) 

(Lyth 1984) providing critical constraints on inflationary models (see Baumann et al. 2008, for a review). It is therefore no 
surprise that a number of sensitive ground-based/sub-orbital CMB polarization experiments are either planned or in progress. 



* Gravitational lensing of CMB B-modes by intervening matter generates a B-mode anisotropy (see Lewis and Challinor, 2006, for a 
review). This effect will be ignored in this paper, since we will concentrate on the detectability of tensor modes at low multipoles I 20, 
where the effects of lensing are small. 
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Examples include BICEP (Yoon et ai, 2006), CLOVER (North et al, 2008), EBEX (Oxley et ai, 2004), QUIET (Seiffert e< 
al, 2006) and SPIDER (Crill a/., 2008). In addition, groups in Europe and the US have considered designs for a B-mode 
optimised space satellite capable of probing tensor-scalar ratios r £ 10~ 2 (de Bernardis et al. 2008; Bok et al. 2008). 

The Planck satellite* is scheduled for launch in April 2009 and has polarization sensitivity in 7 channels over the frequency 
range 30-353 GHz. As described in SPP05, the sensitivity of Planck limits B-mode recovery to low multipoles I < 20*. 
Nevertheless, until a new polarization-optimised satellite is flown, Planck is the only experiment capable of probing these low 
multipoles. It is therefore important to analyse Planck's performance for B-mode detection in the presence of realistic noise 
levels and polarized foregrounds and therefore to assess whether it can provide useful complementary data to experiments 
probing higher multipoles. This is the main goal of this paper. 

The problem of detecting primordial B-modes at low multipoles is unusually difficult. Unlike the detection of temperature 
anisotropies, Galactic foregrounds are expected to have a much larger amplitude than any putative primordial B-mode signal 
over the entire sky (see Section 2). Accurate foreground removal is therefore essential for B-mode detection at low multipoles. 
There has been an enormous amount of work on CMB foreground subtraction (see e.g. the review by Delabroille and Cardoso 
2007; Leach et al. 2008). A variety of methods have been developed, ranging from 'blind' techniques that make few physical 
assumptions concerning the Galactic foregrounds (examples include Internal Linear Combination (ILC), e.g. Bennett et al. 
(2003), Independent Component Analysis e.g. Hyvaxinen (1999) and its fast and spectral- matching variants (Maino et al. 
2002; Delabrouille, Cardoso and Patanchon 2003)), 'semi-blind' methods that make limited use of prior information on the 
foregrounds (such as Maximum Entropy, Stolyarov et al. (2002) and template matching (Bennett et al. 2003; Slosar, Seljak 
and Makarov 2004; Slosar and Seljak 2004; Eriksen et al. 2004a)) and parametric fitting techniques based on physical models 
of the foregrounds (Eriksen et al. 2006; Eriksen et al. 2008; Dunkley et al. 2008a). A number of other methods have been 
developed which incorporate some aspects of these techniques and use, for example, wavelet or harmonic decompositions 
(e.g. Tegmark, de Oliviera-Costa and Hamilton 2003; Hansen et al. 2006; Norgaard-Nielsen and Jorgensen 2008). Some of 
these methods provide approximations to the likelihood function for cosmological parameters (eg. Eriksen et al. 2008) and 
there has also been some work (Gratton 2008) directly addressing the question of modelling the likelihood function from 
multi-frequency maps. Almost all of these methods have been developed for temperature foreground subtraction. In contrast 
to the temperature data, the noise level of Planck polarization maps will be high. There is therefore limited information in 
polarization and hence a restriction on the complexity of polarized foreground removal algorithms that the data can support. 

The problem of B-mode detection in the presence of Galactic foregrounds has been considered by a number of authors. 
Tucci et al. (2005) performed a Fisher matrix analysis for idealised experiments including foregrounds. Amblard, Cooray 
and Kaplinghat (2007) investigated harmonic ILC subtraction for B-mode detection at high mulipoles (£ £ 20) for various 
experimental configurations. Betoulle et al. (2009) considered the application of Spectral Matching Independent Component 
Analysis (SMICA) to perform a Fisher matrix analysis for various experiments, including Planck. The work most closely related 
to ours is the paper by Dunkley et al. (2008b), which focuses on B-mode detection with a future satellite with high signal- 
to-noise in polarization, rather than the low signal-to-noise case relevant to Planck. Low signal-to-noise introduces additional 
complexity to the foreground subtraction problem, nevertheless, there are strong similarities between our approaches. 

The layout of this paper is as follows. The Planck Sky Model, which is used in this paper to model polarized Galactic 
foregrounds, is described briefly in Section 2. This model is compared to realizations of the primordial CMB polarization 
signature to determine the magnitude of the foreground subtraction problem. ILC foreground subtraction is described in 
Section 3. We show that the ILC method is fundamentally limited for B-mode detection because of the offset associated with 
the dominant B-mode signal. Section 4 analyses foreground template subtraction techniques and we present a heuristic model 
for constructing a pixel-based polarization likelihood function. This is applied to simulations with Planck-like noise. Section 
5 introduces a classification scheme for foreground subtraction methods based on their dominant errors. Our conclusions are 
summarized in Section 6. 



2 THE MAGNITUDE OF THE PROBLEM 
2.1 The Planck Sky Model 

The Planck Sky Model (PSM) has been developed by the Planck Component Separation Working Group for use in simulations 
of the Planck mission. Summaries of the model are given by Leach et al. (2008) and Dunkley et al. (2008b) and a detailed 
description will be provided in a forthcoming paper by Delabrouille et al. (2009, in preparation). The polarized foreground 

t For a description of Planck and its science case, see 'The Scientific Programme of Planck' (2005), hereafter SPP05. 
t Although in theory it may be possible to extract some information on primordial B-modes from the noise dominated data at higher 
multipoles, a number of systematic effects such as cross-polar leakage and errors in the polarizer angles are expected to dominate at 
higher multipoles. 
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Table 1: RMS residuals outside internal mask 1 



Map T OK) Q OK) U QiK) 



CMB 


50.4 


0.257 


0.232 


CMB B-modc 


2.26 


0.062 


0.064 


PSM 30 GHz 


97.2 


4.29 


3.94 


PSM 44 GHz 


36.0 


1.51 


1.37 


PSM 70 GHz 


19.0 


0.661 


0.608 


PSM 100 GHz 


22.7 


0.720 


0.692 


PSM 143 GHz 


40.8 


1.29 


1.27 


PSM 217 GHz 


119.2 


3.82 


3.76 


PSM 353 GHz 


874.1 


28.3 


27.9 



1 For maps smoothed with a Gaussian of FWHM 7°. 



model used in this paper is similar to that used by Dunkley et al. (2008b). Briefly, the model includes polarization from a 
power-law synchrotron component with geometrical suppression factors, polarization angles and polarization fractions based 
on the magnetic field model of Miville-Deschenes et al. (2008). It also includes a power-law dust component based on an IRAS 
dust template derived by Finkbeiner, Davis and Schlegel (1999). Polarization from point sources is ignored. The resulting 
dust polarization fraction in this model is ~ 5% over most of the sky, corresponding to the 'high' polarization fraction used in 
Dunkley et al. (2008b). Preliminary indications from the BICEP experiment suggest a significantly lower polarization fraction 
of ~ 1-2% close to the Galactic plane, though this figure may be unrepresentative of regions at higher Galactic latitude where 
depolarization may be lower (BICEP collaboration, private communication). 

Figure 1 shows Q and U simulations at Healpix (Gorski et al, 2005) resolution NSIDE=2048 for a single realization of 
the concordance cold dark matter model§. The model has a tensor-scalar ratio r = 0.1. The PSM of polarized foregrounds 
is shown in the lower panels over the frequency range 70-217 GHz with the colour (grey) scales adjusted to span the range 
—30 fjK to 30 fiK. The structure of the foreground is fairly similar at each of these frequencies, and so we use the 217 GHz 
maps to define a polarization mask by simply applying a threshold to each of the Q and U maps. An 'internal mask' is then 
defined as the union of the two Q and U masks. For the tests described in this paper we use a fairly conservative mask which 
removes 37% of the sky. 



2.2 Impact of foregrounds on B-mode detection 

Figure 1 is not particularly useful for assessing tensor mode detection with Planck, since Planck will be noise dominated in the 
B-mode for all but the lowest multipoles. Smoothed maps, as shown in Figure 2, are of more relevance. Here we show the maps 
of Figure 1 in the regions outside the internal mask after smoothing with a Gaussian of FWHM 7° and repixelization to a 
resolution of NSIDE=16. One can see that at this resolution (almost signal dominated for Planck) the peak-to-peak variations 
in the polarization maps of the primary CMB are of order ~ 0.7 /iK. The PSM is shown for the regions outside the internal 
mask, but unlike Figure 1, the temperature scales of the colour tables are set by the true maximum and minimum values of 
the maps. Evidently at the most sensitive Planck channel (143 GHz) foregrounds dominate over the primary CMB signal over 
the whole sky. In fact, to get a better feel for how accurately we need to subtract foregrounds for .B-mode detection, Figure 
3 shows the contributions to the Q and U maps from the B-mode alone for r = 0.1. The rms contribution from the B-mode 
at this resolution is about a quarter of the rms of the B-mode, and so foreground subtraction to significantly better than 5% 
accuracy at 143 GHz is required to detect a B-mode with r — 0.1. As we will show, this presents a formidable challenge for 
Planck, even with the simplified foreground model assumed in the PSM. The rms values for the PSM and for the primary 
CMB maps shown in Figures 2 and 3 are listed in Table 1. Power-spectra for the maps shown in Figure 2, computed using 
the pseudo-CV (PCL) estimator^, are plotted in Figure 4. 



§ Throughout this paper, apart from the tensor-scalar ratio r, we use the cosmological parameters derived from the 3-year WMAP data 
(Spcrgel et al. 2007) assuming zero curvature and a single scalar spectral index. The tensor spectral index is fixed at Tit = 1. 
" Note that as described in Efstathiou (2006) these PCL estimates are highly sub-optimal at low multipoles, but they are perfectly 
adequate for illustrating the magnitude of the foreground problem. We will discuss more optimal methods in Sections 3 and 4. 
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Figure 1. Q (left) and U (right) maps: the upper panel shows CMB simulations. The remaining panels show the PSM at 70, 100, 143 
and 217 GHz. The temperature scale (thermodynamic temperature) is in fiK. All maps were generated at Healpix NSIDE=2048. 
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Figure 2. As Figure 1, but with maps generated at Healpix NSIDE=16 and a smoothing of 7° FWHM. The internal mask described in 
the text has been applied to the PSM. 
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Figure 3. Maps of the S-mode contributions to the primary CMB realizations (r = 0.1) shown in Figure 2. 




5 10 15 20 25 30 5 10 15 20 25 30 



Figure 4. PCL E and B-mode power spectrum estimates computed for the CMB simulations and foreground components of Figure 
1. The power spectra are computed for the region of the sky outside the internal mask. No instrumental noise has been added to the 
simulations. The blue points show the power spectrum estimates for the CMB. The red lines show the theoretical input CMB spectra. 
The foreground power spectra are as follows: 70 GHz (dark blue); 100 GHz (purple); 143 GHz (orange); 217 GHz (green). 



3 ILC COMPONENT SEPARATION 
3.1 Summary 

The internal linear combination method is very simple and well known (see e.g. Bennett et al. 2003; Eriksen et al. 2004b). 
Suppose we have M maps T l {p) (temparature, Q, or U polarization at each pixel p) at different frequencies, we find the linear 
combination 

M 

T(p) = ^ Wl T» (2) 
that minimises the variance 

Var(T) = (T 2 ) - (T) 2 (3) 
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subject to the constraint that the sum over the Wi is equal to unity. (Angular brackets in this Section denote averages over 
map pixels.) Equation (3) can be minimised using a Lagrange multiplier. If the map covariance matrix 

Ma = ^ - (r))(T\p) - <T J » (4) 

p 



is invertible then the solution is 

T, Mr 1 



(5) 



and we denote the corresponding map T ILC (p). (If Mij is not invertible then a family of solutions for Wi exists.) As a slight 
variation of this method for noisy data, one might consider subtracting the noise variance from equation (4) so that we 
minimise primarily on the foreground residuals. In this case, the solution is identical to (5) with Mij replaced by 

Mij — > Mij — WjSij, (6) 

where <r 2 is the noise variance of the map at frequency i. (Note that this variant is inconsistent with a Bayesian formulation 
of ILC presented in Gratton (2008).) 



3.2 ILC with variable foregrounds 

Let us ignore instrumental noise, but assume that the foregrounds differ with frequency. At each frequency i, we can write 
the data vector as 

T l = S + F i (7) 

where S is the frequency independent CMB signal and F l is the foreground. (For simplicity, we assume that S and F % have 
zero mean). To find the ILC solution, we must extremize the quantity 

S 2 = WiWj((S + F l )(S + F j )} + A j ^ w i ~ 1 J ( 8 ) 

ij \ i / 

with respect to the weights Wi and A. The solution is 

(i + E^^E*^ 1 



V F- 1 

t—tkj kj 

where 



Xi = (SF l ) 



.} 



(10) 



and we have assumed in (9) that there is enough variation in the foregrounds with frequency that the matrix Fij is non-singular 
and therefore invertible. Notice that in the limit of zero signal, Xi — 0, the solution is just equation (5) with Mij = Fij. In 
this case the variance of the ILC solution is 

<(AT ILC ) 2 >= (iX 1 ) • (11) 

This gives an indication of the best possible foreground subtraction achievable through ILC, i.e. of the limitation imposed 
by 'foreground mismatch'. (We will apply the term 'foreground mismatch' generally to subtraction techniques to denote the 
residual contamination by foregrounds even in the absence of CMB signal and instrument noise.) The terms in equation (9) 
proportional to Xi lead to an offset in the ILC solution with an amplitude proportional to Xi and independent of the amplitude 
of the foreground. Such offsets exist even if there is not enough variation in the foregrounds to mimic the CMB. Consider the 
well known case of two frequency channels with an identical foreground template: 



Ti = S + aF, 
T 2 =S + F, 

(without loss of generality, since there are only two unknown functions, S and F) . The ILC solution in this case gives 



(12) 



T lhC = S - g^F, (13) 
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Figure 5. Distributions of the tensor-scalar ratio r for 10 simulations with ILC cleaning as described in the text. The simulations were 
generated with r = 0. 

independent of the parameter a. The error in (13), which we will term 'cross-correlation offset' (called 'cosmic covariance' 
by Chiang, Naselsky and Coles, 2008), is independent of the amplitude of the foreground template and is irremovable, since 
the amplitude of the effect depends on the specific realization of the CMB signal. Equation (10) generalises this result to an 
arbitrary number of channels and includes the bias caused by foreground mismatch. The cross-correlation offset averages to 
zero over a large number of realizations of the CMB. However, it can be significant for any single realization. Cross-correlation 
offset is particularly serious for B-mode polarization measurements. This is because we are trying to detect a B-mode signal 
that is a small fraction of the B-mode contribution to the Q and U maps (c/. Figure 3). The cross-correlation terms Xi for 
the Q and U maps are therefore dominated by the contribution from B-modes, and since these offset terms are fixed for a 
particular realization of the B-modes, they set a fundamental limit on the amplitude of a B-mode signal detectable via ILC 
subtraction. 

This is illustrated by Figure 5, which shows the likelihood distributions of the tensor-scalar ratio, r, for ten ILC cleaned 
CMB realizations generated with r = 0. In each case, we generated an NSIDE=16 low resolution CMB map smoothed with a 
Gaussian of FWHM 7° and added the PSM foregrounds at the four frequencies 70, 100, 143 and 217 GHz as plotted in Figure 
2. The internal mask of Figure 2 was applied to each map and an ILC cleaned map was computed from the solution (5). 

The pixel likelihood function is 



In the examples shown in Figure 5, we compute the likelihood function for the Q and U maps allowing only the parameter r 
to vary. (This is a very good approximation to a full likelihood analysis because r is weakly correlated with other cosmological 
parameters.) A small diagonal noise of 0.1 /J.K was added to the Q and U maps to regularise the signal covariance matrix, 
which is otherwise numerically singular because the maps are over-pixelized. Figure 5 shows that the cross-correlation offset 
causes biases r ~ 0.03, but with a large spread so that biases in excess of r ~ 0.10 are seen. The results shown in Figure 5 
are insensitive to the number of channels used in the ILC. 

The effects of the cross-correlation offset on the power spectra are illustrated in Figure 6 for the specific CMB realization 
of Figure 2 with r = 0. Here we plot the quadratic maximum likelihood (QML) estimates (Tegmark and de Oliveira-Costa 
2001; Efstathiou 2006) 




(14) 



where x is the (T, Q, U) ILC map and C is the sum of the signal (S) and noise (N) covariance matrices 



C i:j = {xiXj), C = S + N. 



(15) 



(16) 



where 



y\ = XiXjElj , 



r=(T,X,E,B), 



(17) 
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Figure 6. QML estimates of the E and B-mode polarization spectra for the realization plotted in Figure 3 with r = 0. The upper panel 
shows the analysis of the CMB maps on the cut sky. The lower panels show the analysis of the ILC cleaned CMB maps. The error bars 
show the diagonal components of (20) using the theoretical input values of r for each realization. The solid lines show the theoretical 
input spectra for r = 0.1. 



with matrices E r given by 

e " = \°- , %k c - 1 - 1181 



x and C are the data vector and data covariance matrix defined in equation (14), and T is the Fisher matrix 

(19) 



dc g _ 1 dc c _ 1 



The covariance matrix of the QML estimates is given by the inverse of the Fisher matrix (19) 

(ACtACfi) =T~ X . (20) 

The upper panels of Figure 6 show the QML estimates for the CMB maps alone computed from the cut sky. (Regularizing 
noise of 0.1 fiK was added to the maps as in the likelihood analysis of Figure 5.) The error bars on the points are computed 
from the diagonal components of the Fisher matrix (19) assuming Cf = 0. The lower panels show the QML power spectra 
for the ILC cleaned maps, illustrating the cross-correlation offset. 

A significant cross-correlation offset is inherent in any 'blind' component separation technique. For the ILC method as 
presented above, the size of the offset is set by the amplitude of the B-mode signal with a resulting large effect on the 
B-mode. In principle, the cross-correlation offset could be reduced by first disentangling E and B-modes and performing 
ILC independently on maps constructed from the two sets of modes. If the whole sky is available, the separation into E and 
B-modes is exact and unambiguous. If the sky coverage is incomplete, one could decompose the maps into almost pure E and 
B-modes following the techniques described by Lewis (2003). Such a decomposition will necessarily lead to information loss 
since ambiguous modes must be discarded. Even if near pure B-modes are identified, so reducing the cross-correlation offset, 
we will see in the next subsection that an ILC solution can amplify the instrument noise to unacceptably high levels. As we 
will show in Section 4, rather than invoke some type of modal decomposition, it is possible to reduce the cross-correlation 
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Table 2: ILC weights 





70 GHz 


100 GHz 


143 GHz 


217 GHz 








Galaxy no noise 


W\ 


U! 2 


W3 


U14 




_ jCMBj2^1/2 


(fiK) 


T 


-0.1427 


-1.0457 


2.995 


-0.8070 




0.947 




Q 


-1.5525 


5.2210 


-2.936 


0.2675 




0.0054 




u 


-1.6387 


5.5929 


-3.285 


0.3303 




0.0047 




Galaxy+CMB 


Wl 


U! 2 


W3 


U14 




_ jCMBj2jl/2 


( M K) 


T 


37.191 


-129.68 


113.52 


-20.028 




9.547 




Q 


-5.3553 


19.526 


-15.707 


2.536 




0.0578 




u 


-6.7704 


26.042 


-21.886 


3.614 




0.0619 





Notes: The upper table lists the ILC weights, tu t , computed from the foregrounds alone (i.e. equation (5) with Fij replacing Mij). The 
column labelled ((T ILC - T CMB ) 2 > 1 /2 u sts t h e rms residual of the ILC cleaned maps and the true CMB maps for the regions outside the 
internal mask shown in Figure 2. The lower table lists the ILC weights for the combined noise-free foreground and CMB maps shown in 
Figure 2. 

offset to negligible levels by template fitting and to avoid catastrophic amplification of instrument noise. The scheme outlined 
in Section 4 has the added advantage of leading to a simple model for the full polarization likelihood function. 



3.3 ILC with noise 

So far, we have ignored instrumental noise. In the presence of instrumental noise, the solution for the ILC weights is approx- 
imately (9) with Fij replaced by 

Fi^Fij+OiSij, (21) 

where af is the noise variance of the map at frequency i. (This solution is exact if the noise-signal (NS) and noise-foreground 
(NF) covariances can be neglected.) If noise dominates in the frequency covariance matrix, the ILC solution corresponds to 
inverse noise-variance weighting, 

— ■ (22) 

which has, of course, nothing to do with foregrounds. One can get a rough idea of how noise will affect the ILC solution by 
computing 

<(AT ILC ) 2 >= (XX 1 ) ( 23 ) 
instead of (10). In the limit that noise dominates, the ILC residuals will simply reflect the noise level of the ILC map, 

fe^ , (24) 



-2 

0"ILC 



where we have assumed that the noise at each frequency is independent. 

Table 2 lists the ILC weights for the foregrounds alone (equation 5 with Fij replacing Mij ) computed for the pixels outside 
the internal mask shown in Figure 2. The ILC residuals using these weights are listed in the last column of the table. For the 
Q and U maps the residuals are an order of magnitude lower than the rms contribution expected for a B-mode with r — 0.1. 
Thus, if the PSM model is correct, it is in principle possible to assign weights that subtract foregrounds to achieve a limit of 
r ~ 10~ 3 (i.e. for the PSM, 'foreground mismatch' is negligible). However, three of the four weights are considerably larger 
than unity and will therefore amplify any instrument noise that is present (equation 24). If somehow we were given these 
weights, we would only able to make use of them if the instrument noise at each channel were extremely low. For example, 
using the instrumental sensitivities in SPP05 for Planck, these weights would lead to disastrously high noise levels of ~ 0.6 
/iK for the ILC cleaned Q and U maps. Such a high noise level is of no use for B-mode detection. 

The situation become even worse if we include the CMB in the ILC solution (lower entries in Table 2). The ILC solution 
listed in Table 2 is for the CMB realization of Figure 2 with r = 0. The polarization power spectra for this solution are shown 
in the lower panel of Figure 6. The cross-correlation offset amplifies the ILC residuals by more than an order of magnitude 
compared to the foreground-only solution (see the last column in Table 2) and so some foreground B-mode leaks into the 
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power spectrum at low multipoles as shown in Figure 6. This is consistent with the fundamental limit of r ~ 0.1 imposed by 
the cross-correlation offset discussed in the previous section. However, notice that the weights are now even larger than for 
the foreground- only case and so any instrumental noise will be highly amplified in the ILC solution. 

For realistic experiments, we are therefore caught between a rock and a hard place. In the presence of instrumental 
noise, we would like to minimise the noise when combining frequency channels (equation 22). However, this will not remove 
foregrounds. If we use weights that minimise the foreground residuals (which we cannot find in principle because of cross- 
correlation offset) we amplify the noise to unacceptable high levels. 



4 TEMPLATE FITTING 

The discussion of the previous Section shows that a purely blind component separation method such as ILC is fundamentally 
limited for B-mode detection by the cross-correlation offset, even if foreground mismatch is negligible (i.e. a linear combination 
exists which eliminates the foregrounds to high accuracy). To reduce the cross-correlation offset, a semi- blind technique is 
required which utilises supplementary information on either the foregrounds or the primary CMB signal. In this Section we 
investigate template fitting and show that this provides an acceptable method for B-mode analysis for Planck. 



4.1 Summary 

Let us model the data vector as 
x = s + F/3 + n, 



(25) 



where s is the signal, F is a template matrix, is a vector of unknown parameters and n is the pixel noise vector. For example, 
the data vector x could be a vector of length 2N P consisting of the Q and U maps x = (Q, U), could be a vector of four 
unknown amplitudes (0q, 0u, 0q, 0u) t , and F a 2N P x 4 matrix consisting of two Q and two U foreground template maps 



( f qW 



Fq(N p ) 












FHl) 



Fu(N P ) 
We find by minimising 



*g(l) 



F 2 Q i N v) 







\ 








Fu(l) 



Fu(N p ) J 



(26) 



X 



= (x-F^) T C- 1 (x-F^), 



where C is the covariance matrix (15). The solution is 
= (F T C- 1 F)- 1 (F T C- 1 x). 

The minimum variance estimate of the signal vector, s, is the Wiener-filtered 
s = SC _1 (x-F/9), 



(27) 



(28) 



(29) 



(see e.g. Rybicki and Press 1992). If the data vector is noise-free and contains zero foreground, we see that template matching 
recovers a 'biased' estimate of the signal, 



s = s-F(F T C _1 F)- 1 F T C" 



(30) 



This is the analogue of (13) for template matching (and is identical for a single foreground/template if the covariance matrix 
C is diagonal). Notice that as with equation (13) the offset is independent of the amplitude of the foreground template. Even 
if there is no foreground in our signal, template matching will produce a cross-correlation offset in the recovered signal that 
is independent of the amplitude of the foreground. As with the ILC method, the cross-correlation offset gives a fundamental 
irreducible limit on the amplitude of a B-mode component detectable by template matching. The critical difference with 
the ILC method is that the amplitude of the offset depends on the mismatch between the foreground matrices, Fq(i)Fq(J) 
etc., and the appropriate sub-matrices Cqq etc. of C. The bigger the mismatch, the smaller the cross-correlation offset. The 
method is therefore 'semi-blind' because it uses prior information on the signal+noise covariance matrix to determine the 
vector 0. As we will show below, this prior information reduces the cross-correlation offset by more than an order of magnitude 
compared to the ILC method. 
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Figure 7. Distributions of the tensor-scalar ratio for template subtracted noise-free combined frequency Q and U maps (frequency range 
70—217 GHz), constructed as discussed in the text. Two curves are shown, one using the 30 and 217 GHz PSM maps as templates and 
the other using the 30 and 353 GHz PSM. 

If the data vector is noise-free, but contains a foreground component F', F' 7^ F/3 the minimum variance signal estimate 
will contain a foreground-mismatch term in addition to the cross-correlation offset 



Evidently, both terms must be shown to be small for B-mode detection. 

As an illustration of the method, we have applied it to the noise- free maps illustrated in Figure 2. We first construct 
inverse noise variance weighted Q and U maps (equation 22) of the CMB (r = 0) + PSM for the four frequency channels 70, 
100, 143 and 217 GHz using the SPP05 detector sensitivies. (The motivation for this will be made clear in the next Section). 
The resulting data vector x, although noise-free, is contaminated with both synchrotron and dust polarization foregrounds. 
We then use the 30 GHz PSM to define a low-frequency foreground template and we use either the 217 GHz or 353 GHz PSM 
to define a high-frequency foreground template. We then compute the template subtracted data vector, x — F/?, and feed this 
into the likelihood function (14) to compute the posterior distribution of r. The results are shown in Figure 7. The 30 GHz 
and 217 GHz PSM templates remove foregrounds to extremely high accuracy of r ^ 10 -3 . The cross-correlation offset term 
in equation (31) is negligible. The primary source of error is caused by the foreground mismatch in equation (31), though for 
these choices of template, the error is small. Using the 30 GHz and 353 GHz templates, the peak of the likelihood is offset by 
~ 3 x 10~ 3 as a result of foreground mismatch term at high frequencies. However, the offset is much smaller than the width of 
the likelihood distribution expected for Planck sensitivities and is therefore ignorable. As we will see in Section 4.3, because of 
Planck's high detector noise levels, it is better to use the 353 GHz channel as a high-frequency polarization template because 
it has a higher foreground signal-to-noise than the 217 GHz channel. 

4.2 Parameter estimation with template marginalization 

In this subsection, we will discuss the problem of parameter estimation from multiple noisy maps with foreground subtraction. 
The aim is to use the results of this subsection to construct an approximate, but easily calculable estimate of the likelihood 
function. 

The data (/, Q and U) per pixel i at each frequency is arranged into a single column vector X. Assume the instrument 
noise is described by a corresponding 'large' covariance matrix N. The foregrounds are modelled by a a mean map at each 
frequency with some spread. We denote the mean by the vector Fo and encode the spread via a Gaussian distribution with a 
covariance matrix P. (P can be very general, being a 'large' matrix, but it can encode 'simple' uncertainties such as a global 
spectral index uncertainty, or an uncertainty about whether a particular pixel is heavily contaminated or not.) 

We assume that CMB has a blackbody, spectrum, so that the CMB-induced signal from each pixel is given by es(i), or, 
for the entire data-set, by Es, where s is the entire CMB signal, and E = I ® e. The matrix E therefore duplicates the CMB 
signal s into each of the frequency bands. Finally, we assume that the foregrounds and CMB are statistically independent. 

Since the CMB is assumed to be Gaussian, the CMB signal s is distributed as 



s = 



s - F(F T C" 1 F)" 1 F T C" 1 s + [F' - F(F T C' 1 F)" 1 (F T C- 1 F')]. 



(31) 
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' e-^ 3 ls ds (32) 



S| 

where S is the signal covariance matrix. The model is assumed to affect the CMB only through S, so the probability for a 
model may be equated to the probability for S. From our assumptions above, the prior on the foreground signal is taken to 
be 

1 e -l(F-F )-P-i(F-F ) dF (33) 



(Note we assume that P 1 exists, though the derivation given below can easily be modified if it does not.) Noise introduces 
a mismatch between X and Es + F, so the probability density for the data given the signal and foreground is, 



J_ -^(X-Es-F) r N" 1 (X-Es-F) 



dX. (34) 



To obtain the probability for a model given the data, we need to combine equations (32)-(34), multiply by a prior p(S) and 
integrate over the signal and uncertainties in the foreground model, 

p(S|X)0C P (S) / / dFds ^ =c -^(X- E s-F)- N -l(X- Es -F) 1 e -I(F-F )-p-l(F-F ) 1 ^s^s ( g 5) 



. |N| V|P| VlSI 

From now on we will drop S-independent factors without comment. Writing the exponent as 

_ 2 y N-E ){ F J +(S F \ N-x + P-Fo- )' (36) 

performing the integral over s and F and using the formulae for block matrices and determinants in Press et al. (1992) the 
result simplifies to the compact expression 

P(S) b -1 Y t (s+m) -i y 

where we have defined M via 

M _1 = E T (N + P) _1 E (38) 
and Y is 

Y = ME T (N + P)" 1 (x-F ). (39) 

Note that equation (38) has a simple interpretation. It is just the usual likelihood function (14) applied to a map Y constructed 
from the data with a noise covariance matrix M give by equation (38). Furthermore, the map (39) is basically an inverse-noise 
weighted map over the frequency channels (but note the 'noise' here is the sum of the detector noise N and the foreground 
uncertainty P and that the best-guess foreground Fo is subtracted from the data at the start). 

This derivation provides the justification for the following simplified model for a polarization likelihood for Planck. In 
Section 4.1, we demonstrated that by using two templates, foreground-mismatch could be reduced to levels that are negligible 
compared to the Planck instrumental noise. It is therefore a good approximation to neglect P in comparison to N. The map 

Y can then be approximated by subtracting the best fit templates (as computed in Section 4.1) from an inverse noise-variance 
weighted set of maps. This justifies the procedure used to construct Figure 7. 

Gratton (2008) presents a conceptually similar technique for constructing a likelihood function from noisy data in the 
presence of foregrounds. Given assumed priors, the scheme effectively averages over linearly weighted combinations of frequency 
channels. The advantage of that technique is that it self-consistently marginalizes over foregrounds, rather than selecting 
specific channels as templates as described above. In fact, in the noise-free limit, to the extent that (25) (with n = 0) is 
an accurate model of the data, the technique gives best-fit weights, u>i, that also yield (30) as the most likely CMB map, 
mitigating foreground mismatch. However, with noise, the algorithm for the resulting likelihood function is considerably more 
complicated. A detailed comparison between the two approaches is currently in progress. 



4.3 Application to simulations with noise 

In this subsection, we apply the above formalism to simulations that include the Planck noise levels. The detector noise 
is assumed to be white and uncorrelated in each pixel. It would be straightforward to generalise this analysis to include a 
more realistic scanning strategy and correlated 'destriping' noise (all that is required is an appropriate model for the noise 
covariance matrices Nij that includes destriping errors). However, correlated errors are expected to be small for Planck and 
so we ignore them in this analysis (for a detailed discussion see Efstathiou 2007). We therefore add uncorrelated white noise 
to the NSIDE — 2048 primary CMB + foreground maps at each of the polarization-sensitive Planck frequencies. The maps 
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Figure 8. Q and U maps: Upper panel shows the noise-free CMB simulations with r = 0.1 for regions outside the internal mask. The 
second panel shows the foreground subtracted noisy reconstructions computed as described in the text. The third panel shows the B- 
mode contribution to the noise-free CMB maps. The lowest panel shows the noisy foreground subtracted reconstruction of this B-mode 
contribution. 

are smoothed with a Gaussian of 7° FWHM and reconstructed at NSIDE — 16. The low resolution noise covariance matrices 
are computed using the small-angle approximation (see Appendix A), which gives an excellent approximation to the true 
covariance matrices. 

In reality, the templates will contain a primary CMB signal 1 1 and so the likelihood approximation is slightly more 
complicated than implied in Section 4.2. As described above, the data vector x, is constructed as an inverse noise variance 
weighted sum over a set of frequency channels. The covariance matrix of this vector is written as 

(xiXj) = Sij + $ij + N i:i , (40) 
II It is unlikely that 'external' templates would be of any value for S-mode analysis. 
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Figure 9. Distributions of the tensor-scalar ratio for two simulations with r = and r = 0.1. The dotted (red) lines show the distributions 
for noise-free and foreground-free simulations for regions outside the internal mask. The solid (blue) lines show the distributions for 
foreground subtracted noisy simulations, as described in the text. 

where Sij, and Nij are respectively the primordial CMB, residual foreground and noise covariance matrices. Now construct 
the data vector 

Yi = xi-Ftfi, ($ = $Q,u),Xi = (Q,U)), (41) 

where the superscript denotes frequency. If the template subtraction removes the foregrounds, the average of (41) over noise- 
realizations is 

<y 4 > = *(!-£>*), (42) 

k 

and if the coefficients /? are independent of the signal the covariance matrix (YiYj) is 

k k 

The solution for /? is found by iteratively minimising (27) with C replaced by (YiYj) and ignoring any weak correlation 
between the solution and the signal. The final data vector Y and its covariance matrix (43) are then used to compute the 
likelihood function (15). The parameters /? are well constrained by the data and so it is a good approximation to keep them 
fixed at their central values. The main contribution of the /3 to the error budget is via the noise term (43). 

In the simulations described here, we construct the data vector x from the four frequency channels 70, 100, 143, 217 
GHz, since there is little additional sensitivity to primary CMB signal in the other channels. We use the 30 GHz and 353 GHz 
channels as templates. The internal mask described in Section 2 is applied to all channels. The resulting noisy foreground 
subtracted maps are shown in Figure 8. The upper panel of this figure shows the noise-free CMB simulations for r = 0.1 for 
the regions that lie outside the internal mask. The second panel shows the reconstruction after foreground subtraction from 
the noisy maps following the procedure described above. There is clearly a very good correspondence between the two sets of 
maps. The third panel in Figure 8 shows the noise-free contribution of the B-mode to the Q and U maps. The foreground 
subtracted reconstruction is shown in the lowest panel. Again, there is a good correspondence between the maps, but the 
reconstructed maps are very noisy. In fact, instrument noise dominates over foreground mismatch. A substantial component 
of the noise comes from the templates because the 30 and 353 GHz channels of Planck are significantly noisier than the main 
'CMB' channels at ~ 100 GHz. 

The likelihood functions for r are shown in Figure 9. The dotted (red) lines show the likelihood functions applied to the 
noise-free CMB maps (though with diagonal 'regularizing' noise applied, as described in Section 3) for the two simulations 
with r — and r = 0.1 for the regions outside the internal mask. These likelihoods are close to the 'best' that could be 
achieved from a low resolution experiment in the absence of foreground contamination. The results from our noisy foreground 
subtracted simulations are shown by the solid (blue) lines. The distribution for the model with r — is peaked close to 
r = 0, so clearly residual foreground mismatch is unimportant. The increased widths of the blue curves are caused by residual 
instrument noise, including the noise in the template channels. All of the results described in this paper assume a nominal 
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Figure 10. QML estimates of the E and B-mode polarization spectra for the two noisy foreground subtracted simulations used to 
compute the likelihoods of Figure 9. The upper panel shows the power spectra for the realization with r = and the lower panel shows 
the power spectra for the realization with r = 0.1. The error bars show the diagonal components of (20) using the theoretical input 
values of r for each realization. The solid lines show the theoretical input spectra for r = 0.1. 



mission lifetime of 14 months for Planck. This ensures that every detector on Planck covers the sky twice (see SPP05). We 
show in a separate paper (Efstathiou and Gratton, 2009) that an extended mission lifetime for Planck to four full sky surveys 
leads to a significant improvement in B-mode sensitivity. 

Finally, Figure 10 shows the QML power spectra (equations 16-20) for the foreground subtracted noisy realizations with 
r = and r — 0.1. As discussed in Efstathiou (2006), the QML estimator eliminates mixing of E and B-modes at low 
multipoles almost perfectly on a cut-sky, and so the B-mode spectrum for the r = realization is indeed close to zero at 
low multipoles. However, this Figure shows clearly that for the Planck noise levels, there is little information in the B-mode 
spectrum at multipoles I 10. 



5 COMMENTS ON FOREGROUND REMOVAL TECHNIQUES 

As summarized in the Introduction, a large number of diffuse foreground subtraction techniques have been discussed in the 
literature. Some of these are designed to recover 'physical' foregrounds, e.g. separating free- free from synchrotron emission. 
Other techniques are designed to tackle the problem described in this paper, i. e. the recovery of the primordial CMB anisotropy 
from foreground contaminated maps. The remarks in this Section apply to this latter class of techniques. 

In the preceeding Sections, we identified two distinct forms of error, which we termed 'cross-correlation offset' and 'fore- 
ground mismatch'. These two types of error provides an intuitively useful way of classifying component separation methods. 
Table 3 summarizes our proposed classification scheme. Any foreground subtraction technique can be placed somewhere 
between the blind and unblind rows of this Table. 

Briefly, a purely blind technique such as ILC can, in theory**, reduce foreground mismatch to negligible levels provided 
there are enough frequency bands to describe the foregrounds. (This is the case, for example, for the weights listed in the first 
three rows of Table 2). However, in a purely blind technique, there is no external information to distinguish between CMB 



. e. a set of weights exists. 
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Table 3: Classification scheme of foreground removal techniques 



Scheme 


cross-correlation offset 


foreground mismatch 


Blind (e.g. ILC) 


Significant 


Small (given enough frequency bands) 


Semi-blind (e.g. template fitting) 


Small 


Small (given enough templates) 


Unblind (e.g. model fitting) 


Small (if model is correct) 


Small (if model is correct) 



and foreground components with similar structure on the sky. The result is a cross-correlation offset that is independent of the 
amplitude of the foregrounds. Any purely blind technique {e.g. harmonic ILC) will show a cross-correlation offset. Unless one 
can isolate pure B- modes on a cut sky (see e.g. Lewis 2003), cross-correlation between the CMB E- modes and the foregrounds 
can produce a potentially serious cross-correlation offset. 

The amplitude of the cross-correlation offset can be reduced if some additional information is provided. We have classified 
template matching as a 'semi-blind' technique because it makes use of some prior information, though it is not based on a 
physical model of the foregrounds. The method requires a model for the signal (primordial CMB) covariance matrix and 
the templates provide a model for the angular distribution of the foregrounds (spectral index variations can be taken into 
account by adding more templates) . As shown in the previous Section it is possible to reduce both cross-correlation offset and 
foreground mismatch to negligible levels. 

A third class of technique attempts to model the foregrounds by fitting a parametric physical model (Brandt et al. 1994; 
Eriksen et al. 2006; Dunkley et al. 2008b). If the physical model is a correct representation of the truth, it is possible to reduce 
both the cross-correlation offset and foreground mismatch to negligible levels. However, this type of technique is limited, in 
practice, by the number of frequency channels available. The number of independent parameters describing the model must 
be less than or equal to the number of frequency channels. For Planck polarization, this limits the number of independent 
free parameters to be < 7 (if the Q and U models are treated independently). This limits the complexity of the physical 
model, limiting the scope for redundancy checks. Of course, if the model is incorrect the method will be limited by foreground 
mismatch. 

It is also useful to consider how foreground subtraction techniques are affected by instrumental noise. Instrumental noise 
in a purely blind technique is, in a sense, 'uncontrollable'. For example, the 'ideal' weights listed in the first three rows of Table 
2 remove foregrounds to high precision. However, if they were applied to noisy data they would amplify the instrumental noise 
to high levels (because many of the weights exceed unity). For Planck polarization, the resulting noise amplification would 
be unacceptable. We have shown in Section 4.2 that in the template matching approach, instrumental noise is 'controllable' 
provided the templates have high signal-to-noise. Instrumental noise is a major problem for model fitting techniques. As far as 
we are aware, nobody has yet developed a model fitting technique that incorporates prior information on the angular variation 
of the spectral indices of the diffuse foregrounds (which vary slowly over the sky) . As a proxy, model fitting is usually done 
by independently fitting parameters in a very coarsely pixelized map. This reduces the effects of instrument noise on the 
estimated parameters, but even then the effects of noise can be limiting. For example, in the analysis of the WMAP 5-year 
polarization data at a resolution of NSIDE = 8 the synchrotron spectral index was computed at a resolution of NSIDE — 2 
(48 pixels over the whole sky) in Dunkley et al. (2008a). Fairly strong (though not unreasonable) priors were imposed to find 
convergent solutions (e.g. the dust spectral index, which is poorly constrained by the data, was kept fixed). Nevertheless, 
Dunkley et al.'s results for the B-mode power spectrum at low multipoles compare well with those from the template cleaned 
maps of Gold et al. (2008). It remains to be seen whether model fitting can perform well for the more difficult problem of 
B-mode recovery for Planck. We hope to report on this in a future paper. 



6 CONCLUSIONS 

In this paper, we have used the Planck Sky Model to assess the impact of foregrounds on B-mode detection by Planck at 
low multipoles. We have analyzed the internal linear combination technique and shown that the offset caused by B-mode 
polarization pattern (cross-correlation offset) leads to a fundamental limit of r ~ 0.1 for the tensor-scalar ratio even in the 
absence of instrumental noise. This is comparable to the sensitivity limit of Planck if foregrounds are neglected. For realistic 
Planck instrument noise, ILC amplifies the noise of the 'cleaned' polarisation maps to unacceptably high levels. Our results 
show that 'blind' techniques such as ILC are unsuitable for detecting primordial B-modes from a future low-noise 'CMBpoP 
mission. 

We have analysed template fitting, using internal templates constructed from the Planck data and devised a scheme to 
approximate the likelihood function (14) from mult i- frequency maps. We have shown that this scheme works well for Planck 
and offers a feasible way of recovering primordial B-modes from dominant foreground contamination even in the presence of 
noise. According to the results shown in Figure 9, Planck, after the nominal mission lifetime of 14 months, could set a useful 



© 0000 RAS, MNRAS 000, 000-000 



18 Efstathiou, Gratton and Pad 



upper limit of r 0.05 if there is no primordial tensor mode and may even detect a tensor mode if r ~ 0.1. This provides a 
useful complement to ground based/sub-orbital experiments which cannot probe these low multipoles (I £ 10). These limits 
probe an interesting part of parameter space (see Efstathiou and Chongchitnan, 2006, for a review). For inflation with a power 
law potential, <p a , the scalar spectral index and tensor-scalar ratio are approximately 

2 + a 4a . a 

n ^ 1 -^T' r "lV' »■«■»•« 8(1 -n.)^, (44) 

where N is the number of inflationary e-folds between the time that CMB scales crossed the 'horizon' and the time that 
inflation ends. There are indications from WMAP and CMB experiments probing higher multipoles for a small tilttt of the 
spectral index n s ~ 0.97 (Komatsu et al. 2008; Reichart et al. 2008). If this tilt is correct, then the last of these equations 
suggests r ~ 0.1 for any a of order unity. For example, for JV « 60 (Liddle and Leach 2003), the quadratic potential (still 
allowed by the data) gives n s w 0.97, r « 0.13, within the parameter range accessible to Planck. Failure to detect a B-mode at 
r > 10~ 2 would put pressure on 'high-field' (0 of order Planck scale) inflation models, in which there has been recent renewed 
interest (Silverstein and Westphal 2008). 

The simplicity of the PSM is a source of concern. Following WMAP, there is quite a lot of information available on the 
polarized synchrotron emission. The low frequency channels on Planck will provide additional information at ,$ 70 GHz. Thus, 
there is considerable scope for redundancy checks at low frequencies, for example, by varying templates and by parameter 
fitting. The dust contribution to polarization is much more uncertain. Neither the level of polarization, nor dust spectral index 
variations are well constrained by current data. We will almost certainly have to wait until Planck flies to assess whether 
polarized dust emission poses a serious problem for B-mode analysis. Planck is heavily reliant on the 353 GHz channel to model 
dust polarization, because it has the highest signal-to-noise on the foreground. There is some limited scope for redundancy 
checks using the 217 GHz channel. However, if polarized dust emission is complex, it may not be possible to achieve a limit 
of r ~ 0.1 with Planck using template fitting or any other foreground removal technique. 
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APPENDIX A: PIXEL NOISE COVARIANCES OF DEGRADED RESOLUTION MAPS IN THE 
SMALL ANGLE LIMIT 

Al Temperature 

Let Xi denote the pixel value in the high resolution map and Xi denote the pixel value in the low resolution map. The harmonic 
coefficients computed from the high resolution map is 

aim = ^2 XiCliYl m (9i), (Al) 

i 

where Q.i is the solid angle of a high resolution map pixel. So, the pixel values in the degraded map are 

X, = aernYem(0i)fe(9s), (A2) 

where fi(9 s ) is the smoothing function applied to the high resolution map. In terms of the pixel values of the high resolution 
map, 

X t = x p {l P Y; m (6 p )Y em (d t )f e (6 s ). (A3) 

imp 

Using the addition theorem for spherical harmonics, the pixel noise covariance of the low resolution map is 

(XiXj) = V (x p x q ) + ^ (— A + ^ n p n q P tl (cos e ip )Pe 2 (cos jq )f tl fe 2 ■ (A4) 

^ — ' 47T 47T 

This expression is time consuming to evaluate, but it simplifies significantly if the noise is diagonal (x p x q ) = a p S pq and if we 
assume small angles. In this case, for a Gaussian smoothing function, equation (A4) simplifies to 

- * E^«p (-ft) ^ (-ft) • ( A5 ) 
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A2 Polarization 

In the case of diagonal pixel noise at high resolution, and for small angles, the polarization covariance matrices (QiQ p ) for 
the degraded resolution maps can be approximated by equation (A5). It is interesting to see why this is so. We will consider 
degraded resolution Q maps (the analysis is identical for U maps). The equivalent to (A3) in obvious notation is 

Qi = \ [ 2Y * W 2Y ^ + - 2Y * { ti -2Y(i)\ + iu p [ 2 y» 2 Y(i) - _ 2 F*(p) _ 2 Y(i)])np/<- (A6) 

ILm p 

Now the addition theorem for the tensorial harmonics is 

i — , /Of 4_ 1 \ !/ 2 

si y;- m (ei,0i) S2 y f - m (e 2 ,0 2 ) = (-l) sl - sl Y es2 (f3,a)e* sl \ (A7) 

m 

where we use the Euler angle conventions of Varshalovich, Moskalev and Khersonskii (1988). Applying the addition theorem 
to (A6) gives 

Qi = \ (^r) V2 [ -2 y «c" 7 + 2 ^- 2 c- 127 ] + iup [ - 2 Y l2 J 2 -< - 2 y,_ 2 e-^] ) fip//. (A8) 

f P 

Now we can write 

2 Y lm = 2V2N e AT(G+ ~G-)e lm4, (A9) 
_ 2 y £m = 2V2N e AT(G+ +G-)e lm4, (A10) 

(see Kamionkowski et al. 1997) and so 

Q< = X) 2 ("i^ - ) 2v / 2iVM?(G + + G")n p ^[g p cos(2a + 27)- Mp sin(2a + 27)]. (All) 

ip 

In the limit I — > co, the prefactor in (All), tends to 

(^i)"W- (^i)i 

and 

GL(5)-^ 4 (Jo( S )±J 4 ( S )), (A13) 
and so (All) becomes 

Qi ~ (^r 1 ) Jo + ^ sin 13 1 2 ^ ' 9p cos ^ 2a + 2t ^ ~ Up sin ^ 2a + 27 ^ • ( A14 ) 

ip 

Summing over £, 

Q 1 f Q 2 \ 

Q% ~ ^^fecos(2a + 27)-u p sin(2a + 2 7 )] cxp ( -^|J , (A15) 

P 

for Gaussian smoothing. Thus, if the q p and ti p are uncorrelated between pixels and are uncorrelated with each other, and 
for small angular separations, it is possible to show after some algebra that equation (A15) gives exactly the same covariance 
matrix as the scalar result of equation (A5) 

«ww « ^E«^ (-J§) (-It) • ( A16 ) 

An analogous derivation applies for the covariance matrix (UiUj) (note that (QiUj) « 0). By comparing with numerical 
simulations, we find that the scalar approximation (A16) is an excellent approximation and is perfectly adequate for the 
smoothing scales adopted in this paper. 
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